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LONG  TERM  GOALS 

The  long  term  goal  of  this  research  is  to  significantly  advance  computational  models  for  multi-scale 
flow  physics  in  geometrically  and/or  hydrodynamically  complex  oceanic  and  coastal  ocean 
environments,  through  refinements  in  the  defined  physics,  domain  definition  and  computational  grid 
resolution.  The  particular  focus  is  on  improved  coupling  of  spectral  wave  models  and  circulation 
models  within  the  framework  of  adaptive,  parallel  scalable,  unstructured  computational  grids  for 
ocean,  coastal  ocean,  inlet,  and  inland  waters. 

OBJECTIVES 

The  objective  of  this  project  is  the  dynamic  and  tight  coupling  of  the  SWAN  spectral  wind  wave 
model  and  the  ADCIRC  circulation  model  using  identical,  unstructured,  large  domain, 
computational  grids.  Such  coupling  leads  to  vast  improvements  in  computational  accuracy  and 
efficiency  compared  to  standard  approaches  for  wave-current  coupling.  Implementing  the  wave  and 
circulation  models  on  unstructured  grids  allows  for  wave  energy  and  flow  processes  to  be  optimally 
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and  locally  resolved.  Large  domains  seamlessly  integrate  wave  and  circulation  physics  from  the  deep 
ocean  to  the  shelf  to  inlets  and  floodplains.  The  need  for  grid  nesting  and  the  associated  problems 
with  overlapping  grids,  interpolation  onto  nested  grid  boundaries,  and  two  way  communications 
between  nested  grids  are  all  eliminated.  The  application  of  identical  unstructured  grids  leads  to 
improvements  in  accuracy  since  the  need  to  interpolate  between  heterogeneous  grids  is  eliminated, 
avoiding  the  associated  aliasing  of  inter-code  solutions  and  forcing  functions.  In  addition,  resolution 
requirements  in  the  wave  and  circulation  models  are  often  similar  due  to  the  response-forcing 
function  interaction  of  the  two  models.  Finally  the  use  of  identical  unstructured  grids  provides 
enormous  efficiencies  in  parallel  computing,  since  inter-code  communications  is  intra-core. 

We  also  want  to  continuously  improve  scalability  in  order  to  be  able  to  simulate  larger  and  better 
resolved  grids  on  the  ever  evolving  parallel  computer  systems.  Furthermore  we  want  to  automate  the 
definition  of  the  unstructured  grids  for  the  next  generation  of  wave-circulation  models  within  the 
framework  of  /^-adaptive  discontinuous  Galerkin  (DG)  based  solutions.  These  methods  optimize  the 
application  of  high  resolution  zones  to  correctly  capture  the  wave-circulation  coupling  within  high 
gradient  surface  water  elevation,  current,  and  wave  transformation  zones  without  over-resolving 
adjacent  areas,  effectively  optimizing  the  accuracy  of  a  given  computation  for  a  specified  cost. 

APPROACH 

We  have  focused  on  three  areas  to  accomplish  the  defined  objectives.  The  first  area  is  integrating  the 
SWAN  wave  and  ADCIRC  circulation  model  in  such  a  way  as  to  maximize  the  accuracy  and 
scalability  of  the  combined  models.  The  goal  is  to  remain  scalable  and  efficient  up  to  10,000’s  of 
processors  and  to  achieve  wall  clock  times  that  allow  for  forecasting.  The  SWAN  and  ADCIRC 
model  components  are  applied  to  the  identical,  unstructured  grid;  share  parallel  computing 
infrastructure;  and  run  sequentially  in  time.  Wind  speeds,  water  levels,  currents  and  radiation  stress 
gradients  are  node-based,  and  therefore  can  be  passed  through  memory  or  cache  to  each  model 
component.  Parallel  simulations  based  on  domain  decomposition  utilize  identical  sub-grids,  and  the 
communication  is  highly  localized.  Inter-model  communication  is  intra-core,  while  intra-model 
communication  is  inter-core  but  is  local  and  efficient  because  it  is  solely  on  adjacent  sub-grid  edges. 
This  communications  paradigm  is  summarized  in  Figure  1.  The  resulting  integrated 
SWAN+ADCIRC  system  is  highly  scalable  and  allows  for  localized  increases  in  resolution  without 
the  complexity  or  cost  of  nested  grids  or  global  interpolation  between  heterogeneous  grids. 

Extensive  verification  against  analytical  solutions  has  been  performed  to  demonstrate  consistency 
and  evaluate  convergence.  Hurricane  waves  and  storm  surge  have  been  validated  for  Hurricanes 
Katrina  (2005),  Rita  (2005)  and  Ike  (2008).  Efficiency  for  operational  models  is  demonstrated  for 
large  grids  which  require  less  than  20  minutes  of  wall  clock  time  per  day  simulation. 

The  second  area  is  evaluating  performance  and  scalability  enhancements  for  Continuous  Galerkin 
(CG)  and  Discontinuous  Galerkin  (DG)  implementations  of  ADCIRC.  Detailed  analyses  of  the 
current  operational  CG  version  of  ADCIRC  indicate  that  scalability  of  the  implicit  code  is  limited  by 
the  global  communications  necessary  for  the  matrix  solution.  Analysis  indicates  that  explicit  CG 
implementations  are  not  only  faster  but  have  even  better  scalability.  Furthermore  investigations  of 
the  scalability  of  DG  based  solutions  show  that  they  are  also  better  than  the  CG  based  implicit 
solutions  along  with  improvements  in  accuracy  on  a  per  degree  of  freedom  basis. 
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The  third  focus  area  is  developing  ^-adaptive  schemes  for  DG  solutions  to  the  shallow  water 
equations  (SWE)  in  anticipation  of  coupling  with  a  DG  based  grid  adaptive  versions  of  the  SWAN 
wave  model.  This  will  improve  the  unstructured  grid  capabilities  for  the  ADCIRC  circulation  model 
by  dynamically  improving  resolution  and  ensuring  that  sufficient  grid  resolution  is  provided  for  all 
relevant  flow  scales.  DG  algorithms  are  particularly  well-suited  for  both  propagation  and  advection 
dominated  problems  with  or  without  sharp  gradients  in  the  forcing  function,  bathymetry,  and/or  flow. 
DG  methods  inherently  preserve  mass  perfectly  on  an  elemental  level,  which  make  them  ideal  for 
coupling  flow  and  transport  models.  Example  problems  include  flows  with  strong  eddies  such  as 
those  issuing  from  and  through  inlets;  flows  through  rapidly  varying  bathymetry  such  as  canyons  or 
deep  scour  holes  in  inlets;  and  flows  with  sharp  fronts  such  as  tidal  bores,  ebb  tide  -  waves 
interactions  in  inlets  and  wetting  fronts.  DG  methods  are  conceptually  similar  to  finite  volume 
methods  although  DG  methods  are  readily  implemented  with  higher-order  bases.  DG  methods  are 
also  ideal  since  non-conforming  h  and  p  refinement  and  adaptivity  can  be  implemented.  In  addition, 
DG  methods  are  significantly  more  accurate  on  a  per  degree  of  freedom  basis  than  CG  methods. 

WORK  COMPLETED 

The  SWAN+ADC1RC  code  has  been  completed  and  the  resulting  tightly  coupled  scalable  parallel 
unstructured  grid  wave-current  model  has  been  thoroughly  validated  (Dietrich  et  al.  2009b).  The 
parallel  implementation  of  unstructured  SWAN  uses  ADCIRC ’s  highly  efficient  and  scalable 
parallel  MP1  based  communications  subroutines.  This  synergy  allowed  for  the  very  rapid 
development  of  the  parallel  unstructured  grid  version  of  SWAN  by  Delft  University.  The 
SWAN+ADCIRC  code  has  been  verified  and  validated  and  is  already  being  used  by  numerous 
researchers  around  the  world.  The  resulting  SWAN+ADCIRC  model  will  be  released  as  a  Gnu 
licensed  combined  code  for  use  as  a  community  model. 

The  efficiency  of  the  CG  solutions  to  the  SWE’s  has  been  a  matter  of  intensive  study.  The  scalability 
and  efficiency  of  CG  based  SWE  solutions  can  be  dramatically  improved  through  explicit 
implementations  (Tanaka  et  al.,  2009).  The  improvements  come  from  the  elimination  of  the  MPI  All 
Reduce  calls  used  in  the  Conjugate  Gradient  matrix  solver  which  requires  all  computational  cores  to 
communicate  with  each  other.  It  is  noted  that  explicit  time  stepping  is  not  a  limitation  on  time  step 
which  is  in  fact  controlled  by  wetting  front  propagation  speeds  and  accuracy  concerns.  We  have 
also  completed  a  comparison  study  of  DG  and  CG  methods  that  finds  that  for  comparable  accuracy, 
both  methods  have  very  similar  computational  costs  (Kubatko  et  al.,  2009a). 

Detailed  studies  of  the  time  step  limitation  of  the  two-dimensional  RKDG  SWE  solutions  were 
performed  to  examine  the  stability  properties  of  these  solutions  (Kubatko  et  al.,  2008).  Semi  discrete 
DG  approximations  using  polynomials  spaces  of  degree  p  =  0,  1,2,  and  3  were  considered  and 
discretized  in  time  using  a  number  of  different  strong-stability-preserving  (SSP)  Runge-Kutta  time 
discretization  methods.  Approximate  relations  between  the  two  dimensional  CFL  conditions  were 
derived  and  numerical  results  verify  the  CFL  conditions  that  were  obtained,  and  “optimal”,  in  terms 
of  computational  efficiency. 

Improvements  have  been  made  for  DG  based  implementations  in  the  robustness  and  accuracy  of  the 
RKDG  SWE  based  wetting/drying  algorithm,  including  proving  convergence  for  a  range  of 
problems  (Bunya  et  al.,  2009a).  The  new  wetting  and  drying  treatment  was  verified  through 
comparisons  with  five  problems  with  exact  solutions  and  convergence  rates  were  examined. 


3 


The  DG  SWE  code  p-  (polynomial  order)  and  standard  h-  (grid)  refinement/adaptivity  options  have 
been  further  refined  and  our  progress  is  summarized  by  Kubatko  et  al.  (2009b).  Further 
improvements  in  DG  efficiency  are  being  made  by  implementing  quadrilaterals.  In  addition  to 
having  advantages  from  a  degree  of  freedom  perspective  for  low  order  interpolation,  quadrilaterals 
have  distinct  advantages  from  the  perspective  of  improvements  in  accuracy  for  a  specific  cost 
(Wirasaet  et  al.,  2009). 

Validation  of  the  DG  SWE  code  has  been  performed  by  studying  tides  in  the  western  North  Atlantic 
(Kubatko  et  al.,  2009b)  as  well  as  in  Pacific  Ocean  inlets.  We  have  also  applied  DG  based 
algorithms  to  hurricane  problems  in  Texas.  The  DG  solutions  are  stable,  robust  and  perform  well. 

RESULTS 

The  details  of  the  SWAN+ADCIRC  model  implementation  are  described  by  Dietrich  et  al.,  2009b. 
We  have  verified  that  the  physics  of  wave  induced  coastal  setup  is  correctly  modeled.  Here  we 
briefly  describe  the  validation  of  the  coupled  model  in  computing  hurricane  waves  in  the  Gulf  of 
Mexico  for  Hurricanes  Katrina  and  Ike.  SWAN+ADCIRC  utilize  the  SL15  grid,  Figure  2  and  3, 
which  has  been  validated  for  applications  in  southern  Louisiana  (Bunya  et  al.,  2009b;  Dietrich  et  al., 
2009a).  This  grid  incorporates  local  resolution  down  to  50m,  and  extends  to  the  Gulf  of  Mexico  and 
the  western  North  Atlantic  Ocean,  Figure  4.  It  includes  a  continental  shelf  that  narrows  near  the 
protruding  delta  of  the  Mississippi  River,  sufficient  resolution  of  the  wave -transformation  zones  near 
the  delta  and  over  the  barrier  islands,  and  intricate  representation  of  the  various  natural  and  man¬ 
made  geographic  features  that  collect  and  focus  storm  surge  in  this  region.  The  SL15  grid  contains 
2,409,635  nodes  and  4,721,496  triangular  elements.  Details  of  all  parameter  selections  for  both 
SWAN  and  ADCIRC  are  given  by  Dietrich  et  al.  (2009b). 

Figure  5  depicts  the  computed  significant  wave  heights  at  12-hr  intervals  as  Katrina  enters  the  Gulf, 
generates  waves  throughout  the  majority  of  the  basin,  and  then  makes  landfall  in  southern  Louisiana. 
In  its  early  stages,  Katrina  generated  significant  wave  heights  of  6-9m  in  the  eastern  half  of  the  Gulf. 
However,  as  the  storm  strengthened  on  28  August  2005,  the  significant  wave  heights  increased  to  a 
peak  of  about  22m  at  2200  UTC,  and  waves  of  at  least  3m  were  generated  throughout  most  of  the 
Gulf.  SWAN+ADCIRC  has  been  validated  to  a  variety  of  wave  and  water  level  measurement  data 
sets.  In  deep  water,  the  National  Data  Buoy  Center  (NDBC)  collected  and  analyzed  wave 
measurements  at  12  buoys.  Figure  6-8  compare  measured  significant  heights,  mean  directions  and 
mean  periods  to  computed  values  from  SWAN+ADCIRC  as  well  as  with  results  from  the  WAM 
model.  SWAN  matches  the  measurements  well  at  most  buoys,  and  in  fact  the  results  are  quite 
similar  to  WAM.  Extensive  comparisons  to  all  data  are  presented  by  Dietrich  et  al.  (2009b) 

Figure  9  shows  the  evolution  of  Hurricane  Ike  while  it  tracked  through  the  Gulf  of  Mexico  with 
significant  heights  reaching  about  19m,  and  large  waves  being  generated  throughout  the  entire  Gulf. 
The  SWAN  wave  solution  for  Ike  has  been  compared  to  measured  results  from  NDBC  buoys.  The 
significant  wave  heights,  directions  and  mean  periods,  shown  in  Figure  10  through  12  suggest  that 
SWAN+ADCIRC  is  performing  as  well  for  Ike  as  it  did  for  Katrina  and  Rita. 

SWAN+ADCIRC  was  benchmarked  on  Ranger,  a  Sun  Constellation  Linux  Cluster  with  62,976 
cores,  at  the  Texas  Advanced  Computing  Center  (TACC)  (http://www.tacc.utexas.edu).  Simulations 
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were  performed  with  the  coupled  model  and  again  with  its  individual  components  in  order  to  discern 
coupling  effects  on  simulation  times.  The  models  were  run  on  256  to  5,120  cores.  As  shown 
in  Figure  13,  the  individual  SWAN  and  ADCIRC  models  both  scale  linearly  through  about  1,000- 
1,500  cores,  but  they  diverge  at  higher  numbers  of  computational  cores.  ADCIRC’s  timing  results 
level  off,  because  the  global  communication  associated  with  its  implicit,  conjugate-gradient  solver 
begins  to  dominate  the  simulation  time.  The  highly  localized  solution  procedure  in  SWAN  allows  it 
to  scale  linearly  through  5,000  cores,  enabling  performance  of  less  than  5  min  per  day  of  Katrina 
simulation.  Also  note  the  coupled  model  shows  linear  scaling  to  about  3,000  computational  cores, 
but  then  it  levels  off.  At  this  point,  the  communication  overhead  from  ADCIRC  implicit  CG 
solution  slows  down  the  coupled  model.  However,  the  performance  in  this  range  is  about  20  min  per 
day  of  Katrina  simulation,  which  is  sufficient  for  forecasts  of  large  storms. 

Tanaka  et  al.  (2009)  investigated  the  scalability  of  implicit  and  explicit  versions  of  the  CG 
implementation  of  ADCIRC.  Figure  14  compares  explicit  and  implicit  simulations  of  the  SL15  grid 
(2,351,832  nodes)  and  a  uniformly  refined  by  a  factor  of  four  version  of  the  grid,  SL15-04x 
(9,314,706  nodes).  We  note  that  the  explicit  simulations  are  about  twice  as  fast  as  the  implicit 
implementations.  More  importantly,  the  implicit  code  shows  weak  parallel  scaling  with  a  limit  in 
scalability  that  is  improved  with  larger  grids.  However,  the  explicit  code  shows  much  improved 
scalability  compared  to  the  implicit  code.  Figure  15  and  16  show  that  this  is  largely  related  to  the 
MPI  All  Reduce  in  the  required  global  communications  in  the  matrix  solver  in  the  implicit  code. 

This  then  suggests  that  the  CG  implementation  of  ADCIRC  be  run  explicitly,  resulting  in  further 
improvements  in  scalability  of  the  SWAN+ADCIRC  code.  A  study  of  the  performance  of  the  DG 
implementation  of  ADCIRC  was  also  done  showing  that  the  explicit  DG  code  performs  much  better 
in  terms  of  scaling  than  the  CG  code  (Kubatko  et  al.,  2009a). 

The  development  of  the  DG  based  implementation  in  terms  of  explicit  time  stepping,  wetting-drying, 
and  h-p  adaptivity  and  are  summarized  by  Kubatko  et  al.  (2008),  Bunya  et  al.  (2009a),  and  Kubatko 
et  al.  (2009b).  Wirasaet  et  al.  (2009)  describe  the  development  and  use  of  quadrilaterals  within  the 
DG  framework  in  order  to  gain  further  efficiencies  in  accuracy  versus  cost. 

Detailed  studies  of  tidal  elevations  and  currents  in  Willapa  Bay,  WA  show  that  the  DG  based 
solutions  perform  well  in  areas  characterized  by  its  large  tidal  range,  fast  currents,  strong  advection, 
and  large  tidal  flats.  We  have  also  applied  DG  based  solutions  to  compute  storm  surge  during 
Hurricane  Ike  along  the  Texas  coastline,  Figure  17  and  18.  The  computations  are  stable  and  robust. 

IMPACT/APPLICATIONS 

This  work  will  significantly  improve  the  accuracy  of  the  computed  physics  for  coupled  coastal  ocean 
wave  and  circulation  computations.  The  physics  will  be  directly  improved  through  the  much  tighter 
wave  -  circulation  model  coupling,  as  well  as  through  the  ability  to  apply  much  higher  resolution 
grids  which  will  in  fact  refine  themselves  to  the  appropriate  levels  though  the  DG  adaptivity.  This 
work  will  also  reduce  model  development  times  since  the  base  grids  developed  by  hand  only  need  to 
represent  the  geometry  and  not  the  hydrodynamics.  The  grids  will  adjust  themselves  automatically 
and  optimally  to  account  for  the  waves  and  currents  as  the  hydrodynamics  evolves  in  the  simulation. 
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NEIGHBOR  CORE 


Figure  1.  Schematic  of  parallel  communication  between  models  and  cores.  Dashed  lines  indicate 
communication  for  all  nodes  within  a  sub-grid,  and  are  inter-model  and  intra-core.  Solid  lines 
indicates  communication  for  the  edge-layer-based  nodes  between  sub-grids,  and  are  intra-model 

and  inter-core. 
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Figure  2:  ADCIRC  SL15  model  domain  with  bathymetry  (m). 
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Figure  3:  ADCIRC  SL15  bathymetry  and  topography  (m),  relative  to 
NAVD88  (2004.65),  for  southern  Louisiana. 
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Figure  4:  ADCIRC  SL15  grid  resolution  (m)  in  southern  Louisiana. 
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Figure  5:  Hurricane  Katrina  significant  wave  height  contours  (m)  and  wind  speed  vectors  (m  s'1) 
at  12-hr  intervals  in  the  Gulf  of  Mexico.  The  six  panels  correspond  to  the  following  times:  (a) 
2200  UTC  26  August  2005,  (b)  1000  UTC  27  August  2005,  (c)  2200  UTC  27  August  2005,  (d) 
1000  UTC  28  August  2005,  (e)  2200  UTC  28  August  2005  and  (f)  1000  UTC  29  August  2005. 
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MEAN  WAVE  DIRECTION 


DATE  IN  2005 


Figure  7:  Mean  wave  directions  (°),  measured  clockwise  from  geographic  north,  during 
Hurricane  Katrina  at  12  NDBC  buoys.  The  measured  data  is  shown  with  black  dots,  the  modeled 
SWAN  results  are  shown  with  black  lines,  and  the  modeled  WAM  results  are 

shown  with  gray  lines. 
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SOLUTE  MEAN  WAVE  PERIOD!; 


DATE  IN  2005 

Figure  8:  Mean  wave  periods  (s)  during  Hurricane  Katrina  at  12  NDBC  buoys.  The  measured 
data  is  shown  with  black  dots,  the  modeled  SWAN  results  are  shown  with  black  lines,  and  the 

modeled  WAM  results  are  shown  with  gray  lines. 
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Figure  9:  Hurricane  Ike  significant  wave  height  contours  (m)  and  wind  speed  vectors  (m  s'1)  at 
12-hr  intervals  in  the  Gulf  of  Mexico.  The  six  panels  correspond  to  the  following  times:  (a)  1800 
UTC 10  September  2008,  (b)  0600  UTC 11  September  2008,  (c)  1800  UTC 11  September  2008,  (d) 
0600  UTC  12  September  2008,  (e)  1800  UTC  12  September  2008  and 
(f)  0600  UTC  13  September  2008. 
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SIGNIFICANT  WAVE  HEIGHTS  (m) 
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Figure  10:  Significant  wave  heights  (m)  during  Hurricane  Ike  at  12  NDBC  buoys.  The  measured 
data  is  shown  with  black  dots,  and  the  modeled  SWAN  results  are  shown  with  black  lines. 
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Ml vAN  WAVE  DIRECTIONS  (°) 


Figure  11:  Mean  wave  directions  (°),  measured  clockwise  from  geographic  north,  during 
Hurricane  Ike  at  12  NDBC  buoys.  The  measured  data  is  shown  with  black  dots,  and  the  modeled 

SWAN  results  are  shown  with  black  lines. 
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ABSOLUTE  MEAN  WAVE  PERIODS  (s) 
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Figure  12:  Mean  wave  periods  (s)  during  Hurricane  Ike  at  12  NDBC  buoys.  The  measured  data 
is  shown  with  black  dots,  and  the  modeled  SWAN  results  are  shown  with  black  lines. 
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Figure  13:  Timing  results  for  SWAN+ADCIRC  and  its  components  on  the  TACC  Ranger 
machine.  The  times  shown  are  wall-clock  minutes  per  day  of  Katrina  simulation  on  the  SL15 
grid.  SWAN  results  are  shown  in  red,  ADCIRC  results  are  shown  in  blue,  and  SWAN+ADCIRC 

results  are  shown  in  purple. 
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Figure  14:  Performance  as  a  function  of  the  number  of  cores  of  the  SL15  grid  implicit  (with  1  sec 
time  step),  SL  15x04  grid  implicit  (with  0.5  sec  time  step),  SL15  grid  explicit  ( with  1  second  time 
step),  and  SL15x04  grid  explicit  (with  0.5  second  time  step)  ADCIRC-CG  simulations  on  Ranger. 
Wall  clock  time  is  in  seconds  for  a  0.25  days  tidal  simulation. 
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Figure  15:  Component  performance  of  the  SL15-04x  grid  with  implicit  ADCIRC-CG  simulation 
on  Ranger.  Wall  clock  times  are  in  seconds  per  one  day  of  simulation  time  and  normalized  to  1 
second  time  steps.  Components  examined  include  total  time;  the  continuity  (GWCE)  equation; 
momentum  equation;  wet/dry  evaluation;  local  communication  between  overlapping  layers  of 
adjacent  subdomains;  global  communication  used  the  MPI  all  reduce  associated  with  the  implicit 

code  matrix  solution. 
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Figure  16:  Component  performance  of  the  SL 15 -04x  grid  with  explicit  ADCIRC-CG  simulation 
on  Ranger.  Wall  clock  times  are  in  seconds  per  one  day  of  simulation  time  and  are  normalized  to 
one  second  time  steps.  Components  examined  include  total  time;  the  continuity  (GWCE) 
equation;  momentum  equation;  wet/dry  evaluation;  local  communication  between  overlapping 

layers  of  adjacent  subdomains. 


22 


Figure  1 7:  ADCIRC  DG  computed  water  column  height  (m)  during 
Hurricane  Ike  along  the  Texas  coastline. 
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Figure  18:  ADCIRC  DG  computed  velocity  (m/s)  during  Hurricane  Ike  along  the  Texas  coastline. 
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